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ABSTRACT 

Lightning is a major cause of damage in laminated composite aerospace 
structures during flight. Due to the dielectric nature of Carbon fiber reinforced 
polymers (CFRPs), the high energy induced by lightning strike transforms into 
extreme, localized surface temperature accompanied with a high-pressure 
shockwave resulting in extensive damage. It is crucial to develop a numerical tool 
capable of predicting the damage induced from a lightning strike to supplement 
extremely expensive lightning experiments. Delamination is one of the most 
significant failure modes resulting from a lightning strike. It can be extended well 
beyond the visible damage zone, and requires sophisticated techniques and 
equipment to detect. A popular technique used to model delamination is the 
cohesive zone approach. Since the loading induced from a lightning strike event is 
assumed to consist of extreme localized heating, the cohesive zone formulation 
should additionally account for temperature effects. However, the sensitivity to this 
dependency remains unknown. Therefore, the major focus point of this work is to 
investigate the importance of this dependency via defining various temperature 
dependency profiles for the cohesive zone properties, and analyzing the 
corresponding delamination area. Thus, a detailed numerical model consisting of 
multidirectional composite plies with temperature-dependent cohesive elements in 
between is subjected to lightning (excessive amount of heat and pressure) and 
delamination/damage expansion is studied under specified conditions. 
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INTRODUCTION 


Carbon fiber reinforced polymers (CFRPs) offer significant weight and 
performance advantages over metals, making them ideal for aerospace applications. 
However, CFRPs are dielectric (compared to highly conductive metals, such as 
aluminum); thus when subjected to an electrical impulse (such as a lightning strike) 
the electricity is not conducted through the CFRP structure. Consequently, the 
energy from the lightning strike is transformed to thermal energy through resistive 
super-heating of a highly localized region of the structure near the arc attachment 
point(s). The extreme surface temperatures induced by the lightning, along with an 
accompanying high-pressure shockwave (caused by rapid super-heating of the air 
surrounding the arc) yields extensive visible damage near the arc attachment point, 
as well as damage away from the arc attachment point that may not be detected 
through visual inspection. The strike zone (red zone, Figure 1) is the major visible 
damage zone after a composite panel is hit by a lightning strike. Excessive surface 
temperatures cause major visible damage such as fiber sublimation and matrix 
pyrolysis (material recession) and delamination in this area. This strike zone and 
some area beneath this zone (the blue repair zone in Figure 1) are fully removed 
and replaced during regular inspection procedures. However, a major area of 
concern is the residual zone (yellow zone, Figure 1), which is not repaired or 
replaced after inspection, but might contain highly-damaged areas due to lightning. 
As a result of a lightning strike, some regions in the residual zone are exposed to 
very high temperatures (much higher than matrix Tg), causing major degradation of 
matrix, fiber, and the interface. Furthermore, exposure of the matrix to such 
temperatures may increase the susceptibility of the composite to creep deformation 
during normal service hours. Consequently, delaminations and inelastic matrix 
deformations are among the major damage mechanism occurring in this zone. 
However, these sub-critical damages are not easily detected during regular visual 
inspections, thus they may not be fully included in the repair zone. Accumulation of 
damage in overlapping residual zones, over multiple lightning strike events, might 
unexpectedly impact the structural performance of the composite panel 
significantly. Based on numerical predictions of the extent of the accumulated 
damage in the residual zone, the repair zone (blue zone, Figure 1) might have to be 
redefined, so that so that repairs can be completed properly, avoiding any 
unexpected catastrophic failure. 



Figure 1. Damage resulting from simulated lightning strike experiment (courtesy of NASA 
Langley Research Center (LaRC)) 


Although a predictive numerical tool can provide valuable information on the 
extent of the lightning-induced delamination damage in the residual zone, 
developing a competent numerical tool is a very challenging task. Available studies 
regarding simulated lightning strikes on composite panels are extremely limited 
[1,6], with almost no attempt to model the damage induced by lightning in CFRPs. 
Meanwhile, the variation of CFRP properties with high temperatures induced by 
lightning and the influence of these properties on lightning-induced damage 
remains largely unknown. Starting with one of the most extensive available 
experimental studies, a series of simulated lightning strike experiments was 
conducted by Feraboli and Miller [2] on un-notched and filled-hole carbon/epoxy 
laminates. The damage area was thoroughly evaluated via ultrasonic scanning, and 
advanced optical microscopy to gain further understanding of the morphology of 
damage [2], Later Kawakami and Feraboli [3] evaluated the effect of repair 
procedure on the structural performance of wire-mesh protected scarf-repaired 
carbon/epoxy specimens following a lightning strike. According to the authors [3], 
a poor repair performs equally or worse than a fully unprotected specimen, and 
therefore the damage margin (the region to be repaired) has to be defined precisely. 
Yet, another experimental study by Hirano et al [4] has addressed the evolution of 
damage in graphite/epoxy composite laminates due to lightning strikes with 
variations in lightning parameters and specimen size. Ultrasonic testing, sectional 
observation and micro X-ray inspections were used to assess the damage, which 
was mainly categorized into fiber damage, resin deterioration, and internal 
delamination modes. Their results showed that the lightning tests created large 
delaminations propagating in the shape of a pair of fans starting from the lightning 
attachment point in each interlayer. In addition to delamination, matrix cracks were 
also observed, however delamination was emphasized as the dominant failure mode 
among the others. The experiments in [4] also demonstrated that each damage mode 
shows correlations with a particular lightning (loading) parameter, nevertheless, 
variation of specimen geometry has little effect on the damage response. Later 
Ogasawara et al [5] carried out a finite element (FE) coupled thermal-electrical 



analysis of carbon fiber reinforced polymer composites (CFRP) exposed to 
simulated lightning current, in order to investigate the damage behavior caused by 
the lightning strike. In their study, the delamination area was estimated to be the 
same as matrix resin decomposition area (the areas where the temperature was > 
300 C), as delamination was assumed to occur because of matrix resin 
decomposition accompanied by dielectric breakdown of pyrolysis gases. However, 
no actual continuum delamination modeling was accomplished within the FEM 
simulation [5], 

Based on the mentioned literature [2-5], delamination is the most significant 
failure mode induced by a lightning strike, which is meanwhile very hard to inspect 
through regular inspections. When inter-laminar damage (delamination) is the 
dominant failure mechanism, one recent appealing technique used in capturing 
delamination is the cohesive zone approach [6-11]. The cohesive zone approach 
uses interface elements incorporating various traction-separation laws to capture 
crack initiation and propagation. Since the lightning phenomenon is assumed to be 
a combination of extreme temperatures and a pressure shockwave leading to 
extensive inter-laminar damage, the cohesive zone formulation should additionally 
account for temperature effects. A comprehensive numerical model to predict the 
amount of lightning-induced delamination in the residual zone is currently not 
available in the mentioned literature. As the physical phenomenon behind a 
lightning strike is a very complicated one, a step-by-step methodology shall be 
followed to devise a general modeling tool for lightning-induced delamination 
prediction. The first step, which is the major focus of this work, is a systematic 
sensitivity study to present the effect of temperature-dependent parameters on the 
delamination progress, utilizing a temperature-dependent cohesive zone 
formulation. 

Following this introduction, constitutive equations of the temperature-dependent 
cohesive user element are explained concisely. Next, construction of the numerical 
model, required input parameters, and various patterns of temperature dependency 
for these parameters are depicted. Numerical simulations of the conducted 
lightning-induced delamination are presented, and the obtained results are discussed 
after. Finally, a brief summary and conclusion is presented. 


NUMERICAL MODEL 

The numerical FE model consists of individual elastic, transversely isotropic 
plies and interface elements in between them as shown schematically in Figure 2. 
As in-ply damage is considered to be less significant based on available 
experimental evidence [1-5], the only damageable portion is considered to be the 
interface. The mathematical damage model used for modeling interfacial damage in 
this work is implemented as a User Element in the FE code ABAQUS by 
Naghipour et al. [12], and is described briefly in this section. 



Figure 2. Schematic of lamina and interface elements in the FE model 

The pursued mixed-mode cohesive element formulation is based on the 
constitutive models suggested by Pandolfi and Ortiz [10] and Turon et al. [11] and 
is further modified to consider the temperature dependency of interface properties. 
Figure 3 shows the global and local geometries of the three-dimensional 8-node 
zero-thickness interface element, which are related by the standard isoparametric 
mapping. Nodes 1-4 represent lower face of the interface and nodes 5-8, which 
coincide geometrically with nodes 1-4, represent the upper surface. (The zero 
thickness has been offset for better visualisation). 
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Figure 3. Interface element with (left) global system and (right) local system 


With respect to the three-dimensional global coordinate system, x, y, z, each 
node has three degrees of freedom, u, v and w respectively. Meanwhile, in the local 
coordinate system the vector of relative displacements between each pair of the 
corresponding upper and lower nodes of the element is defined as: 
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n indicates the opening or normal component and s and t indicate the two shear 
directions, respectively, as shown in Figure 3 (right). The mixed-mode damage 
formulation is based on constitutive interface tractions (in, is, it) and relative 
displacements (5 n , 5 S , 5t) (Figure 4). 
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Figure 4. Mixed-mode traction separation law for the cohesive element 

Interfacial damage initiates upon satisfaction of the quadratic interfacial traction 
interaction criterion as shown in Equation (2). The relative displacement 
corresponding to damage initiation in the mixed-mode plane, 8m° (Error! 
Reference source not found.), is called the mixed mode opening displacement of 
the fictitious crack tip. 
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in 0 , is 0 , it 0 are the normal and shear elastic limits of the interface, which are 
assumed to vary linearly with temperature (0). Once the damage is initiated, it 
starts evolving based on an energy based propagation criterion given by 
Benzeggagh and Kenane (B-K) [13]: 

G,c{e) + {G IIC (9)-Gi c m^fr^r = G c (0) (3) 

G t (0) =Gj(0) +G shear {d ) 

Note, that the form of both the initiation criterion (Equation (2)) and the mixed- 
mode propagation criterion (Equation (3)) are based upon empirical observations. 
It is assumed that these forms are adequate for this material system, and testing this 
hypothesis is beyond the scope of this work. Future experimental studies will be 
conducted to determine if Equations (2) and (3) can sufficiently fit experimental 
data under a variety of temperatures and mixed-mode conditions. G/ c (0) and 
GiiciQ ) are temperature-dependent fracture toughness values for mode I and mode 
II that must be determined from fracture experiments conducted under various 
temperatures. The mixed-mode failure response of the material is described by 
plotting the total critical fracture toughness G c vs. different mode mixities 
( G shear/' G t) . Parameter r| in Equation (3) maintains the shape of the failure locus in 
the mixed mode plane, and the most accurate failure criterion is the one matching 
the material response when plotted in this mixed mode diagram. As the damage 
evolves based on the mentioned energy-based criterion, an isotropic damage 


variable, d, (Equation 4 ) degrades the interfacial tractions, till the final separation 
point ( 8 m f in Figure 4 ) is reached. The damage variable (d) degrading the interfacial 
tractions is calculated at the end of each increment and is compared to determine 
the initialization and the extent of propagated damage in interface elements. 
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FE MODEL DESCRIPTION AND INPUT MATERIAL PROPERTIES 

The numerical model was created using the object oriented ABAQUS Scripting 
Interface (ASI, python) in ABAQUS for further optimization and parametric 
studies. The laminate is a 500x500 mm square composite panel, made of 10 
AS4/PEEK plies (1.4 mm thickness) with a quasi-isotropic stacking sequence and 
interface elements placed in between laminae. Each lamina is defined using elastic, 
transversely-isotropic, 8 node, reduced integration, continuum shell elements. The 
interface elements, implemented as user-defined elements (UEL) in ABAQUS, are 
placed in between the plies to capture the delamination behaviour. 

The material properties for the lamina and the interface (for room temperature) 
are listed in Table I. Lamina properties, longitudinal and transverse ply stiffness 
(En, E22), Poisson’s ratio (v 12), longitudinal and transverse shear modulus (G12, 
G23) and their variation with temperature are partially reported in literature [15], 
Longitudinal and transverse ply thermal conductivities (K11, K22, K33) are taken 
from previous literature [4] and assumed to be constant here. However, in order to 
have a precise numerical model, directional thermal-conductivity data should also 
be determined as a function of temperature through standardized conductivity 
experiments. In contrast to ply properties, nothing is available on variation of 
interfacial properties with temperature in reported literature, Typical interface, 
which presumably vary with temperature, are normal and shear fracture toughness 
(Gic, Giic) and corresponding strength values (x° n , x°s). The fracture toughness values 
at room temperature are measured through available fracture experiments [12], 

In order to thoroughly investigate the effect of temperature dependence of the 
delamination process, the modeling study has been divided into two categories. In 
the first category, three different temperature dependent profiles for interfacial 
fracture/ strength values have been assumed. Single mode and mixed mode fracture 
toughness and strength values (G and x) experience a linear reduction to 80%, 40% 
and 20% of their room temperature values, as the temperature rises form room 
temperature to 4000 °C. In the second category, mixed mode fracture toughness and 
strength values experience a linear reduction to 80%, 40% and 20% of their room 
temperature values, when the temperature rises from room temperature to the glass 
transition temperature (Tg -190 °C) For temperatures higher than Tg, the 
corresponding toughness and strength values are assumed to be close to zero. Note 
that, assuming different temperature dependent profiles with higher/lower rates of 
reduction will undoubtedly yield different delamination responses. The profiles here 
are just initial assumptions, and are not representative of any experimental data. The 
main focus is to emphasize the significance of temperature-dependent interfacial 


properties on lightning-induced delamination response via employing a generic 
temperature dependence profile. In order to determine precisely the temperature- 
dependent fracture toughness values, fracture experiments must be conducted under 
various temperatures, and corresponding fracture loads and crack propagation data 
must be recorded. However, obtaining fracture data in elevated temperature ranges 
is a challenging task. Currently, no experimental standards exist for this, and no 
comprehensive experimental studies have been conducted with definitive success. 
The results reported in a few studies addressing high temperature fracture testing 
[16, 17] seem to be inconsistent and need further investigation. Hence, the results of 
this numerical study can avail as a firm incentive for further research planning on 
temperature-dependent fracture experiments. 


Table I. LAMINA/INTERFACE PARAMETERS AT ROOM TEMPERATURE 
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When the panel is subjected to a lightning strike, the transition from the applied 
electric load into a combination of mechanical (pressure) and thermal (excessive 
heat) loads is a complicated task. As a result of the resistive joule heating caused by 
the dielectric nature of the CFRP, the front surface temperature is estimated to be 
about 3900-4200 °C [4], However, the numerical simulations presented in 
Ogasawara (2010) contained a maximum applied temperature of 3000 °C, which is 
the fiber sublimation temperature. Above these temperature the composite material 
has completely ablated, and cannot transfer any heat; therefore, a maximum applied 
temperature of 3000 °C was also used in the simulations presented in this work. 
Meanwhile a 50 kA strike is supposed to create a high pressure shockwave with an 
approximate 20-30 MPa pressure in the vicinity of the attachment point [1, 14]. In 
order to verify the accuracy of the applied loads, the applied temperature and 
pressure are calibrated so that the backside temperature and displacement correlate 
with the data from available experimental lightning strikes on unprotected panels. 
Applied loads are calibrated successfully to match maximum backside temperature 
(~180 °C) and maximum displacement (~30 mm) of the panel (Figure 5). 
Accordingly, loading boundary conditions (excessive surface temperatures and a 
pressure), are applied to the panel in a circular area (radius=10 mm depth= 0.4 mm) 
in the middle of the panel (Figure 6), and the edges are clamped. The load is applied 
in a quasi-static manner, as for now there is no possibility to extract the damage 
variable from an explicit user-defined element in ABAQUS. However, dynamic 
effects are planned to be included in future studies. As a result of this study number 
of damaged interface elements is recorded and resulting delamination profiles under 
the above-mentioned specified boundary conditions are presented. 




Figure 5. Applied pressure and temperature boundary conditions (calibrated to match the measured 
backside displacement and temperature of the panel) 



Figure 6. Example FE mesh of the composite plate and applied boundary conditions 




RESULTS AND DISCUSSION 


The importance of using accurate temperature-dependent interface properties to 
predict the lightning induced delaminations is demonstrated in Error! Reference 
source not found.. Three various delamination profiles are compared with each 
other; first one with no temperature dependency, and the second one with an 80% 
temperature-dependent linear descent in interface strength and toughness, when the 
temperature rises to 4000C, and the third one with an 80% temperature-dependent 
linear descent in interface strength and toughness, when the temperature rises to Tg. 
While we notice a homogenous propagation scheme in the first case, semi fan- 
shaped delamination propagation is observed for the second case, and we notice a 
totally random damage profile, distributed over a much wider area in the third case. 



Figure 7a. Delamination profiles: No temperature dependency is assumed (left), 80% linear descent 
in interfacial properties (strength and toughness) is assumed (left) 
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Figure 7b. Delamination profiles: No 80% linear descent in interfacial properties (strength and 
toughness) is assumed when the temperature rises to Tg 


The number of damaged elements in the delamination area does not exhibit a 
significant rise in the first and second case (Figure 8), however, the delamination 
profiles are completely different from each other. For the third case, there is a 
drastic increase in the number of damaged elements as well (Figure 8), which 
emphasizes the influence of a fast reduction rate of interfacial properties on 
lightning-induced delamination. Based on observed delamination schemes, the 
mentioned repair zone (blue zone, Error! Reference source not found.) might 
have to be redefined to avoid any structural failures. 
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Figure 9. Comparison of number of damaged elements for various reduction rates of fracture 

toughness in interface elements 
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CONCLUSION 

A numerical model, consisting of multidirectional plies with temperature- 
dependent cohesive interface elements in between, was subjected to lightning 
simulated as a combination of a mechanically applied pressure together with very 
high thermal gradients. A parametric study was carried out, i.e., corresponding 
interfacial parameters were varied over a definite temperature range, in order to 
understand their influence on the predicted lightning damage response. As a basis 
for comparison, the delamination (number of damaged interfacial elements) was 
recorded under the specified conditions. It is concluded that, interfacial fracture 
toughness is the primary parameter controlling lightning-induced delamination 
propagation. The size of the delamination is strongly dependent on the reduction 
rate of fracture toughness values with temperature, as higher reduction rates 
produce a significantly greater number of failed interface elements. It is also 
observed that the shape of the delamination zone is a significant function of 
temperature-dependent fracture properties. Based on the shape of the delamination 
zone, regular repairs might have to be reconsidered to avoid sub-critical failures. 
Therefore, it can be concluded that measuring the temperature dependence of the 
fracture toughness is critical in order to predict the delaminated area resulting from 
a lightning strike. It is worth mentioning that this study is a preliminary study, and 
further complexities will be built up systematically to capture all the necessary 
physics and devise a robust lightning-induced damage prediction model. 
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